arXiv:1504.01935v2 [math.NA] 21 Jan 2016 


Double obstacle phase field approach to an inverse problem 
for a discontinuous diffusion coefficient 

Klaus Deckelnick* Charles M. Elliott^and Vanessa Styles^ 


Abstract 

We propose a double obstacle phase field approach to the recovery of piece-wise 
constant diffusion coefficients for elliptic partial differential equations. The approach 
to this inverse problem is that of optimal control in which we have a quadratic fidelity 
term to which we add a perimeter regularisation weighted by a parameter a. This yields 
a functional which is optimised over a set of diffusion coefficients subject to a state 
equation which is the underlying elliptic PDE. In order to derive a problem which is 
amenable to computation the perimeter functional is relaxed using a gradient energy 
functional together with an obstacle potential in which there is an interface parameter 
e. This phase field approach is justified by proving T—convergence to the functional 
with perimeter regularisation as e ^ 0. The computational approach is based on a finite 
element approximation. This discretisation is shown to converge in an appropriate way 
to the solution of the phase field problem. We derive an iterative method which is 
shown to yield an energy decreasing sequence converging to a discrete critical point. 
The efficacy of the approach is illustrated with numerical experiments. 


1 Introduction 

Many applications lead to mathematical models involving elliptic equations with piece-wise 
constant discontinuous coefficients. Frequently the interfaces across which the coefficients 
jump are completely unknown. A common approach for the identification of these coef¬ 
ficients is to make observations of the field variables solving the equations and use these 
values in an attempt to determine the coefficients by formulating an inverse problem for 
the coefficients. This is generally ill posed and in applications it is usual to use a fidelity to 
the observations functional together with a regularisation of the coefficients. In this paper 
we use a regularisation of the coefficients by employing the perimeter of the jump sets of 
the coefficients. 
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1.1 Model problem 

To fix ideas we consider the following model elliptic problem: 

= 0 in (1-1) 

= ^ on dil, (1.2) 

where is a bounded domain in (d = 2, 3), ^ is given boundary data with zero mean 

[ 9-0 (1.3) 

JdQ 

and a is an isotropic diffusion (conductivity) coefficient. We suppose that the diffusion 
coefficient takes one of the r positive values ai,...,a^. Our interest is in modelling a 
geometrical inverse problem concerning the determination of the regions in which the ma¬ 
terial diffusion coefficient takes these values. Our problem then is to determine the sets 
Ei = {x ^ Q\a{x) — ai} given observations of the solution y of the elliptic boundary value 
problem ( |1.1| ), ( |1.2[ ). In the case of r = 2, under constraints on the nature of the domains 
and boundary conditions, uniqueness and stability results have been proved in [TIE]. In 
this context see also [28] . 

A standard approach is to minimise a fidelity functional 

Jfid{£) \ \yS ~ 2/o6s||o 

over an appropriate class of partitions £ — of where yg denotes the solution of 

the state or forward equation with diffusion coefficient a{x) = G Ei^i = 

l,...,r. Furthermore, O is an appropriate space of observations and yobs G O is given. 
In general this problem is ill-posed and is typically regularised by adding a Tikhonov 
regularisation functional. A numerical approach without regularisation is proposed in [28l 

[32|. 


(aVy) 

dy 


1.2 Geometric regularisation 

In this setting it has been considered appropriate to use perimeter regularisation, [Ml EH 

r 

Jreg{€) ^ H ft), 8 = iEi)U, 

i=l 

where the regularisation parameter a is positive. Minimisers of 

JiS) ■= JfidiS) + Jreg{S) 

are then typically sought in the set of Caccioppoli partitions into r components, i.e. par¬ 
titions £ — of with T-L^{Ei n Ej) = 0,i ^ j, \ IJi=i ^0 ~ ® which 

Ui := XEi belongs to = 1,... ,r. Thus, a Caccioppoli partition corresponds to a 

function u = (i^i,..., Ur) G BV (fl, {ei,..., e^}), where ei,..., are the unit vectors in 
M^. We can then write the regularisation functional in terms of u as follows: 

r 

Jreg{^) — ^ 

1=1 
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Here, \Dui\ is the total variation of the vector-valued Radon measure Dui. Before we 
rewrite the fidelity term let us introduce the Gibbs simplex 

r 

E := {y G |yj > 0,i = 1,... ,r, = 1} 

i=l 

and observe that ei,..., are the corners of S. Consider the set 


X := {u : ^ I u is measurable and u(x) G D a.e. in 

endowed with the L^-norm and define for u G X 


a{u) :='^aiUi (1.4) 

i=l 

and by S{u) the solution of ([Ol) with diffusion coefficient a(u). 

We set a = for later convenience. The constant tt/S arises from the form of the phase 
field relaxation used in ([Tlf , see ( |2.7[ ). 

Problem (PGR) is then to seek minimizers of the functional J : X ^MU{oo} given by 


J(u):= 



if u G BV{yt^ {ei,..., e^}) H X; 


otherwise. 


In this problem the fidelity term is non-convex because of the nonlinearity of the state 
solution operator *S(-) with respect to the coefficient a(u). Also a feature of this natural 
geometric regularisation approach is that the regularisation functional is non-convex. This 
is reflected in the fact that u only takes one of the values which leads to a 

non-convex constraint. 


1.3 Double obstacle phase field approach 

We shall consider a suitable phase field approximation of the above regularisation which 
involves gradient energies and functions that map into the Gibbs simplex. In this approxi¬ 
mation we relax the non-convex constraint u(x) G {ei,. .., by introducing the set 

/C := {u G H^{n,R'^)\u{x) G E a.e. in 

and approximate J by the sequence of functionals Jg : X —> R U {oo}, e > 0 with 

7 / N \ - VobsWo + (^ [ {^\Du\^ + - \uf))dx ,ifuG/C; 

Je(u) 2 Ja^2 2e ( 1 . 5 ) 

I oc , otherwise. 

Here, = J2i=i Iq and for u G X we have /^(l — |u 

Ui). Problem (PDO) is then to seek minimisers of J^. We refer to this approach as a double 
obstacle phase field model because of the constraints 0 < < 1 on the components of the 

phase field vector u. The parameter e is a measure of the thickness of a diffuse interface 
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separating two sets on which the diffusion coefficient is constant. The Cahn-Hilliard type 
energy 



is well established as an approximation of the perimeter functional, see e.g. [121 [m [6]. Note 
that the regularisation remains non-convex through the quadratic Cahn-Hilliard functional 
even though the constraint set is convex. Let us remark that such a phase field model has 
recently been used in a binary recovery problem, see jl5j . 

Note that we view (PGO) as having just one regularisation parameter a. The 6: parameter 
in (PGO) may be viewed as a way of providing an approximation of (PGR) which is 
computationally accessible. 

1.4 Other approaches 

There have been attempts to solve the recovery problem without regularisation of the inter¬ 
faces across which the diffusion coefficients jump. Formally one can write down variations 
of the fidelity functional with respect to variations of the interfaces. For example see [28] . 
In particular the interfaces can be associated with particular level sets of level set functions 
which have to be determined. We refer to [Ml Ea Ea Ei for numerical implementations. 
The use of level set descriptions of the interfaces in the context of perimeter regularisations 
is described in laiMiEii. Related to this is the use of total variation of a regularised Heavi¬ 
side function with argument being a level set function, [2211101. In EH] the authors consider 
the distributed control of linear elliptic systems in which the control variable should only 
take on a finite number of values. To this purpose they introduce a combination of and 
L^-type penalties whose Fenchel conjugates allow the derivation of a primal-dual opti¬ 
mality system with a unique solution. A suitable adaption of this approach could be an 
alternative way to attack the inverse problem considered in the present paper. 

In the different context of image segmentation parametric description of curves have been 
used in conjunction with perimeter regularisation, [8[ l3^. 

On the other hand [iziiasiiM] use total variation regularisation and relax the constraints 
that the indicator functions take just two values. 

1.5 Applications 

Our model problem is an example of the identification of a coefficient in an elliptic equation. 
This problem arises in many applications. For example, a fundamental issue in the use 
of mathematical models of ffow in porous media is that the geological features which 
determine the permeability are unknown. In geology a facies is a body of rock with specific 
characteristics. In our model problem y is the pressure or hydraulic head associated with 
a ffuid (for example, oil or water) occupying the reservoir or acquifer and a is the 
permeability of the rock. We assume that the permeability is isotropic and is piece-wise 
constant. The domains Ei — {x ^ Q. \ a{x) = ai}^ i = 1, 2,..., r model the decomposition of 
the reservoir into facies whose location is unknown. The goal is to use observations of 
the pressure to determine the geometrical decomposition of the reservoir with respect to 
these facies, [2511M1ES]. 
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Such problems also arise in imaging. For example, electric impedance tomography, MM 
US], is the determination of the conductivity distribution in the interior of a domain using 
observations of current and potential. Here y is the electric potential and a is a conductivity 
which takes different values in unknown interior domains. In medical imaging the shape 
and size of interior domains may be inferred from the variation of the conductivity. 

1.6 Outline and contributions of the paper 

• In Section 2 we introduce the functionals and prove that they F-converge to J. 
Furthermore, we show that has a minimum and derive a necessary first order 
condition. This establishes that problems (PGR) and (PDO) have solutions. 

• The optimisation problem in Section 2 is infinite-dimensional. In order to carry out 
numerical calculations we employ a finite element spatial discretisation. This is de¬ 
rived in Section 3 and we prove convergence results for absolute minimizers and crit¬ 
ical points as the mesh size tends to zero. This establishes that the inverse problems 
(PGR) and (PDO) can be approximated by something computable. 

• Section 4 is devoted to formulating an iterative scheme for finding critical points of the 
functional associated with the discrete optimisation problem. The method is based 
on a semi-implicit time discretisation of a parabolic variational inequality which is 
a gradient flow for the energy. In this finite dimensional setting we prove a global 
convergence result for the iteration. 

• Finally in Section 5 we illustrate the applicability of the method with some numerical 
examples. 


2 Problem formulation 


2.1 State equation 

Let Q C be a bounded domain with a Lipschitz boundary. We suppose that g G 
satisfying (1.3) and yobs ^ O are given functions. Here, (O, (•, ‘)o) is a Hilbert space with 
the property that is compactly embedded in O. Furthermore we assume that the 

following Poincare inequality 

||r? - Mom < Cp\\Vrj\\, r? G H\n) (2.1) 

holds, where || • || denotes the L^(fi) norm and Mo{r]) denotes the mean value of 77 with 

Mo{ri) (??, l)o/||l||o, rjeO. 

Typical examples are O — L^(fi) or L?{d^) representing either bulk measurements or 
boundary observations of the solution of the state equation. 

For a given u G X we denote hy y — S{u) G the unique weak solution of the 

Neumann problem 


— V • (a(u)V7/) = 0 in 
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on dCl 


( 2 . 2 ) 

(2.3) 



with M.o{y) — MoiVobs) ill the sense that 

/ a{w)Vy • Vrjdx = / grjdo V77 G 

Jo. JdQ 

Here, a(u) is given by ( |1.4[ ), where we note that 

Cimin ^ ^ Cimax i^ ^5 Uniformly in u G X, 


(2.4) 


(2.5) 


where a^m min(ai,..., a^-), amax •— max(ai,..., a^-). Observe that 5 is a nonlinear 
operator because of the bilinear relation between a(u) and y in ( |2.4[ ). Using ( |2.1| ) together 
with the fact that M.o{y) = Mo (yobs) we infer that the solution y = S{u) satisfies 

Ill/ll < \\y - Mo{y)\\ + M^Mo{yobs)\ < C^W^yW + ^bobsWo. 

\\M\o 


If we combine this estimate with the choice 77 = ^ in (2.4) and use (2.5) as well as the 
continuous embedding ^ L‘^{d^) we deduce that 

ll'S'(u)||iji(a) ^)(ll5'llL2(9a) + llyofesllo) uniformly in u € X (2.6) 

We see that the problem of observing y given u is well formulated because 

S : X ^ O is continuous 

which is a consequence of the following lemma. 

Lemma 2.1. S : X ^ is continuous. 

Proof. Let u G X and {uk)keN ^ sequence in X with ^ u in L^(i4,M’^),fc ^ (X). Since 
0 < Uk^i < l,i = we may assume by passing to a subsequence if necessary that 

u/;, ^ u in and a.e. in Q. Abbreviating y — S{u)^yk = S{uk) we have for 

/ a{uk)X{yk-y)-Xydx= / {a{u) - a{uk))Xy - Xydx. 

Jq Jq 

Choosing y = yk — y we deduce with the help of ( |2.5[ ) and 0 

(^min l|V(yfc-y)|| < ^max |ufc - u\^\Vy\^dx^ ^ 0, fc 


(X) 


by the dominated convergence theorem because 

\uk — up|V 7 /p ^ 0 a.e. in 14, \uk — u\‘^\Vy\‘^ < r\Vy\‘^ a.e in 14 and \Vy\‘^ G L^(14). 
Since Mo{yk — y) — Q we deduce with the help of that S'(ufc) ^ yk ^ y ^ S'(u) in 
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2.2 F—convergence and existence of minimizers 

The use of in the minimization of J is justified by the following F-convergence result. 
Theorem 2.2. The functionals T-converge to J in X. 


Proof. Let us write Je(u) = G(u) + (jFe(u), where G(u) = i||*S(u) — yobsWo continuous 


as a consequence of Lemma 2.1 and the embedding of into O. In Theorem 6.1 

the Appendix we show that 


m 


Fe ^ F, where F(u) = < 



oo 


if u E nX; 


otherwise. 


I—n r 

Using Remark 1.7 in HU we infer that ^ G aF = J. 

Theorem 2.3. The minimization problem minv^x has a solution E 1C. 


(2.7) 


□ 


Proof. Let {uk)keN C /C be a minimizing sequence, Je{^k) \ i^fv^x Since {uk)keN 

is bounded in there exists a subsequence, again denoted by {uk)keN^ sind E 

such that 


Uk Ue in Uk in and a.e. in Ct. 

In particular, E /C. Lemma [2T| implies that S{uk) F(ue) in O which combined with 
the weak lower semicontinuity of the F^-seminorm shows that is a minimum of J^. □ 


Corollary 2.4. Let {u^)e>o be a sequence of minimizers of Then there exists a sequence 
Ck ^ ^ oo and u E BV{PL] {ei,... ^Cr}) HX such that ^ u in L^{Ct^W) and u is 

a minimum of J. 

Proof. By Corollary |6.2| in the Appendix there exists a sequence Ck ^ O^k ^ oo and 
u E BV{Q] {ei,... ,e^}) H X such that ^ u in L^{Q^W). It is well-known that the 
F-convergence of to J implies that u is a minimum of J. □ 


2.3 Necessary first order condition for the phase field recovery 

In order to derive the necessary first order conditions for a minimum of we consider fC as 
a subset of M^). Similarly as in [9], Section 3, one can prove that the solution operator 

S : D 1C ^ is Frechet differentiable with y = F'(u)w,w E 

being given as the solution of 

I a{u)Xy • Vydx = — f a(w)VF(u) • Vrydx Vry E (2.8) 

Jn Jn 


with Mo{y) = 0. As a result, Js is Frechet differentiable on X C L^{Ct^W) n 
with 

J'(u)w = (F(u) - yobs, *5'(u)w) ^ • Fw - -u • w)dx (2.9) 

Jn ^ 
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for w G n In order to avoid the evaluation of *S'(u)w in (2.9) we 

work as usual with a dual problem: Find p G such that Aio{p) = 0 and 


[ a{u)Vp-Vpdx = {S{u) -yobs^v)o V 77 G ( 2 . 10 ) 

JQ 

where we note that the solvability condition {S{u) — Pobs^ l)^ = 0 is satisfied. As a result 
we obtain from (2.9), (2.10) and (2.8) 

J'(u)w = I a{u)\/p •'V[S'{u)w]dx + a f (^eDu • Dw —-u • w)dx 

Jq Jq ^ 

= — a{'w)'VS{u) • Vpdx + a / (eDu-Dw - u-w)dx. 

Jn Jn ^ 


At a minimum u of we have J'(u)(v — u) > 0 for all v G /C. Since a(v — u) = a(v) — a(u) 
we therefore define: 

Definition 2.5. (Phase field critical point) Find u G /C such that for all v G /C 

a I (eDu • D(v — u) — -u • (v — u))dx — [ (a(v) — a(u)) VS'(u) • Vpdx > 0. (2.11) 

Jn ^ Jq 


Remark 2.6. A natural strategy to construct solutions of ( 2 . 11 ) and hence to find can¬ 
didates for at least a local minimum of is to consider the following parabolic obstacle 
problem: Find u(-,t) G /C,t > 0 such that u(-,0) = uq and 

(ut,v —u)+cr I ( 6 :Du • D(v — u) —-u • (v — u))dx — j (a(v) — a(u)) VS'(u) • Vpdx > 0 

Jq ^ Jq 

for all V G /C and all t > 0. Here, p is the solution of (2.10) for u(-,t) and uq G /C is a 
suitably chosen initial function. 

Inserting v = u(-,t — At) into the above relation, dividing by At and sending At ^ 0 we 
formally find that 

||ut||^ + J'(u)ut < 0, 

SO that ^Jg(u(-,t)) < 0 and the value of the objective funtional decreases during the 
evolution. If lim^^oo '^(•5 ^) =: '^00 exists, we expect u^o to be a solution of (2.11). 


3 Finite element approximation 

In what follows we assume that is a polygonal (d=2) or polyhedral (d=3) domain. Let 
us denote by {Th)o<h<ho ^ regular triangulation of Q and set 

Vh^ixe C\n) I X|T e PiiT) for all T e %} C H^Vt) 


as well as 


■■= {x e (VhY I X(^) e E, X G n} c /C. 









Using the construction of the Clement interpolation operator ( 120 ]) it is not difficult to see 
that for every u E /C there exists a sequence (u/i)o</i</io with E fCh such that 

u/i ^ u in as /i ^ 0. (3.1) 

Furthermore, let {yobs)^<h<ho be a sequence of functions E O such that 

Vobs Vobs in e> as 0. (3.2) 


For Uh E JCh we denote by yh = Sh{uh) E Vh the solution of 


(3.3) 


with MoiVh) ^ -^oiVobs)- Here gh : dQ 
tion to g satisfying 


/ a{uh)Vyh ■ Vxdx = / gtxdo Vx G Vh 
Jn Jdfi 

^ is a piecewise linear, continuous approxima- 
/ g^do = 0 and g^ ^ g in L^i^dVt) as /i ^ 0. (3.4) 

JdQ 


In the same way as in (2.6) one can prove that 

II'S'/j(u/,)||hi < c{\\gh\\L^(dn) + WVobsWo) < c uniformly in u/, € ICh, (3.5) 


where the constant c is independent of h in view of (3.2) and (3.4). 

Lemma 3.1. Let {hk)keN be a sequenee with lim^ij^^oo = 0 and E with ^ u 
in L^{n,W). Then Sh^{uh^) S{u) in oo. 

Proof. Let Uk = Uhj^^yk = Sh}.{uk) and y = S{u). By passing to a subsequence if necessary 
we may assume in addition that ^ u a.e. in Lt. Choose a sequence yk E Vhj^ such that 
yk ^ y in Using (2.1) we deduce 

Wvk - Vkllm < Wvk -Vk- Moivk - yk)\\ + M^\Mo{y^bs ~ yk)\ + - yk)\\ 

< c||V(yfc -yfc)|| + {\Mo{y^t - yobs)\ + \Mo{y - yk)\) 

< c||V(yfc - yk)\\ + c{\\y'^^^ - yobsWo + \\y - ykWm)- 

In order to estimate the first term we write 


(3.6) 


/ a{uk)V{yk-yk)-Vxdx 

Jn 

= / a{uk)V{y-yk) ■Vxdx+ / (a(u) - a(ufc)) Vy • Vx^a: + / {ghu - 9)xdo 

Jq Jq JdQ 

for all X E If we let X — Vk ~ Vk ^md take into account (3.6) we obtain 

llyfc-yfcllifi < c\\y-yk\\H^+c(^J^\uk-u\‘^\Vy\‘^dxy 

+c(ll5fefc — 51^2(9^) + \\yoL ~ yobsWo) o,k ^ oo 
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Using Sh we define the following approximation : /C/^ ^ M of J^: 


Je,hi^h) ■= ^\\Shiuh) - VobsWo + <^ + ^(1 - \uhf))dx. (3.7) 

Theorem 3.2. There exists G tCh sueh that Je^hi^h) — Je,h{^h)- Every 

sequenee {uhf^)keN = 0 has a suhsequenee that eonverges strongly in 

and a.e. in Q. to a minimum of J^. 

Proof. Since is finite-dimensional, the existence of a minimum of is straightforward. 
Next, let Uk G /C/,^ be a sequence with lim^^oo hk = 0 and Je^hki^k) = Je.hki^h)- 

Since {uk)keN is bounded in there exists a subsequence, again denoted by 

{uk)keN^ aiid u G /C such that 

Uyjj; ^ u in ^ u in L^{Q^W') and a.e. in f]. (3.8) 

Furthermore, Lemma [3.1| implies that 

Sh.iuk) ^ Siu) in H\n). (3.9) 


We claim that u is a minimum of J^. To see this, let v G /C be arbitrary and irk G JChu a 


sequence with ^ v in see (3.1). Since J^^hki^k) ^ Je.hki^k) we deduce from 

(3.8), (|3.9D and again Lemma [3^ that 


Jg(u) < liminf Je,/j^^(ufc) < limsup Je,fe^^(ufc) < lim Je,hk{^k) = JeiV^ 

k^oo k^oo k^oo 


SO that Je(u) = minvGX: -^e(v). Furthermore, by repeating the above argument with a 
sequence G such that ^ u in we infer in addition that 


lim Je,hk{^k) = ^e(u). (3.10) 

k^oo 


We use this relation to show that llDukW^^ 


|L)u|p. Namely, let us write 


ae 

y 


/ \Dukfdx = Je,ftfc(ufc) [ (1 “ \^k\‘^)dx - i||-S,j,(ufc) - 

^ “ \uf)dx - i||5(u) - VobsWo = y / 1-^ 


= \Du\ dx 


in view of (3.10), (3.8), (3.9) and (3.2). Hence ^ u in H^{n^W) and the theorem is 
proved. □ 


In practice, rather than trying to locate a global minimum of looks for admissible 

points Uh that satisfy the necessary first order condition 


J!,^hi^h){vh - Uh) > 0 for all v/i G/Cfe. 


(3.11) 
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A calculation analogous to (2.11) leads us to the following variational inequality: 

a I {£Duh-D{vh-Uh)--Uh-ivh-Uh))dx- [ (a(v/j) - a(uft)) Vy/j • Vp/idx > 0 (3.12) 
Jn ^ Jn 

for all Vh G where = Sh{uh) and ph G Vh with M.o{Ph) = 0 is the solution of the 
discrete adjoint problem: 


/ 

Jq 


a{uh)Vph ■ Vxdx = {yh - y^bs,x)o Vx € I4. 


(3.13) 


Theorem 3.3. Let {uhj^)ken be a sequenee of solutions of (3.12) with limk^^hk = 0. 
Then there exists a suhsequenee that eonverges strongly in and a.e. in Lt to a 

solution u of (2.11). 

Proof. Let us abbreviate Uk = Uhj^^yk = Shf.{uk) and denote by pk G 14,^ the solution of 


(3.13) with Uh = Uk and yh — yk> Using ( |3.5[ ) and testing (3.13) with x = Pk infer that 
WVkWH^ + WPkWn^ < c uniformly in A; G N. 


Next, inserting ^ (3.12) we deduce 


/ 9 a 

ae / \Duk\ dx < — 

Jq ^ 


[ \uk\^dx+ [ (a(ufc) - i Waj) Vyfc • Vpfcdx 

Ja Jq ^ Gi 


ar , 


< —|fi| + c||Vyfc|| llVpfcll < c. 

Hence, there exists a subsequence, again denoted by {uk)keN^ and u G /C such that 

Uyt ^ u in Uyt ^ u in and a.e. in fl. (3.14) 


Lemma 3.1 implies that 


Vk = Shki^k) S{u) =: y in W{n). 


(3.15) 


Let p G Aio{p) = 0 be the solution of (2.10). Choose pk G 14^ with M.o{Pk) = 0 

such that Pk ^ p in. and write 

/ a{uk)\/{pk - Pk)= / a(ufc)V(p-pfc) • Vydx 

Jq Jq 

+ / («(u) - a{uk))Vp ■ Vxdx + (yk - y, x)o - {y^^s ~ yobs,x)o 
Jq 


for all X G 14^. By choosing X — Pk — Pk and using (2.1), (3.15) and (3.2) we deduce 
\\Pk-Pk\\H^ <4Pk-p\\H^+c( [ \uk-u\'^\Vpfdx)^+c{\\yk-y\\o + \\yobs-yobs\\o) ^0 


which implies that pk ^ p iia H^{Lt). 

Let us next show that u satisfies ( |2.11 ). Given v G /C there exist s a s equence Vk G Khj^ 
such that Vyjj; ^ V in H^{Lt^W) and a.e. in Lt. Then we have from (3.12) 


a 


{eDuk’D{irk-Uk)--Uk’{<rk-Uk))dx- [ {a{irk) - a{uk))Vyk’Vpkdx > 0. (3.16) 
} ^ Jq 
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In order to examine the second term we write 


/ (a(vfc) - a(ufc)) Vyfc • Vpfcda: - / {a{v) - a{u)) V y ■ V pdx (3.17) 

JQ JQ 

= / {a{vk)-ci{^k))[^{yk-y)-^Pk + ^y-^{pk-p)]dx 

JQ 

+ / ((a(vfc) - a(v)) - (a(uA:) - a(u))) \/y • Vpdx ^ 0, ^ oo 

Jn 

since yk ^ y^Pk ^ P in H^{^) where we used ag ain th e dominated convergence theorem 
for the second integral. By passing to the limit in (3.16) and observing that \Du \‘^dx < 
liminf/;;^oo iDukl^^dx we infer that u satisfies (2.11). 


Let us finally show that ^ u in Choose a sequence G fCh^ such that 


^ u in Inserting into (3.12) we obtain 

/ 

Jn 


ae I \Duk\‘^dx < ae I Duk-Dukdx — 


/ Duk'Dukdx-- / Uk-{uk-i^k)dx- (a{uk)-a{uk)) Vyk-Vpkdx 
jQ ^ jQ Jq 


SO that ( |3.14[ ) and ( |3.17[ ) with imply that 

limsup / iDukl^^dx < / \Du\‘^dx. 
k^oo Jq Jq 

Hence j^\D\Xk\^dx ^ In dx, so that Duk —> Du in L^. 


□ 


4 An iterative scheme 

4.1 Iterative method 


Let us consider the following iteration, which can be s een as a time discretization of the 

G ICh let 


parabolic obstacle problem introduced in Remark 2.6l Given G ICh l^t G JCh be 


/( 


the solution of the problem 

- K) ■ (v/i - -Tn [ ia{vh) - a(u))+^)) • Vp^dx (4.1) 

I Jq 

+r„<T / {eDu^'^^ ■ D{vh - - -^h ' (v/i - u)(+^))4x > 0 Vv/, G Kh, 

Jci £ 

where Tn > 0, = Sh{u^) and p'^ G Vh solves the discrete dual problem 

[ ««) • Vx4a; = {y^ - x)o Vx e 14 with MoiPh) = 0- (4-2) 

Jq 

Note that is the unique solution of the convex minimization problem 


mm 1 - 


- <|p - Tn j^a{vh) Vyl ■ Vphdx + TnO jj^^\D\h?dx - • \h)dxY 


12 











4.2 Convergence of the iterative method 


The following result shows that the objective functional decreases in the iteration provided 
the time steps satisfy a suitable condition. In order to formulate it we define 

c := inf I h \ r]e H^{n) \ {0},Mo{r]) = o}- (4-3) 


a : = 


Note that c > > 0 in view of (2.1). 

Lemma 4.1. The sequence satisfies 

n € No, 

provided that 

rn<li + —\\^yl\\LA^Pl\\L^ + ^h\^yl\\l^ 


Or] 


^min 


n G No- 


Proof. Inserting % = into (4.1) we obtain after some calculations 


(4.4) 


-lluy^ - uir + y P(ur - or + y l|uO “ 


< 


+<y + r 1 - + 4(1 - K|"))dx 

/ oW+'jV*; . Vp'idx - f a(ul)Vyj; . Vp'idx =. I + II. (4.5) 

jQx Jn 


Using ( |3.3[ ) for and 2/^^^ with test function p'^ as well as ( |3.13[ ) we may rewrite II as 
follows: 


II = 


[ a«+i)VOi • Vpldx 
Jn 


(4.6) 


= - / a«+i)VO / a«+i)V0^.V(0 

Jq Jq 


= -(o^ - y^ob.s,yr')o + / «(<+rvo^ • v(pr^ - K)dx ^ n,+ih. 

Jq 


Using again (3.13) we may write 

Ih = -jlli/r' - A.\\o + \\\yk- !/lllo - l\K*' - Vtwl - tor' - i.Kjo 

= -^Itor' - A.Ho + - !/lllo - ^Itor' - vlWo 


f 

Jq 


a{ul+^)Vy^ ■ Vpl+^dx, 


,n+l^ 


while 


Ih^ / a«)Vy^V(py^-p)()dx. 

Jq 
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Inserting the above identities into (4.6) and combining it with (4.5) we obtain 


+ J)ll“r‘ - <11" + f PK+i - OIP + l|K+‘ - viwl + JaK+‘) 


< Vfe«)+ / {o,{ul) - a{nl^^))Vyl-V{pl^^-pl)dx 
JQ 

,.n I 


>r'))vs;.‘ ■ vw+‘ 

< j,,iK) + a||V!,j||i»K+‘ - <11 ||VK+‘ - oil. 

It remains to estimate ~ Ph)\\' To begin, note that 

/ aiul+^)Vipl+^-pl)-Vxdx= [ (a«)-a«+l))Vp^Vxdx + (y-+l-y;J,x) 

Jn Jn 

for all X ^ T/i- Inserting x = ~ Ph deduce that 

amm||V(p^+^ -Ph)f 

< d\\Vplh^ 11^+^ -<|| ||V(p;:+^ -Pl)\\ + -PlWo, 


(4.7) 


o 


which implies in view of (4.3) 


iiv(pr‘ -ps)ii < -^iivpsik- ii<+' - <11+^;^iipr‘ -pjiio. 

^min V C ^min 


,,n+l 


Inserting the above bounds into (4.7) and using (4.4) we infer 
1 


u 


"+‘ - <||" + ^||D«+‘ - Oil" + hp;+' - pJllo +- J.MO 


< 


< 


-||Vy;(||Loo||Vp)(||Loo||u”+^ -<f + —^—^l|Vy;(||Loo||y”+^ -y;(||o||u”+^ -< 
in ^min V C 


2 

l|2 


a 1 


,n+l 


n+l 


'mm 
2 


llvpJ||L«||vp;i|i« + h_i||vpf||" 


^min 


LOO 


|<+‘ - u: 


n II2 


Ml +§IK+‘-Pkllo 


< (4-i)iK+‘-<f+hpr‘-pjii?„ 

‘n ^ 


and the result follows. 


□ 


Corollary 4.2. Let G tCh- Then the time steps Tn in (4-1) can he ehosen in sueh a way 
that > 7 > 0,n G N, where 7 depends on the data and possibly on h. For this ehoiee 
the sequenee generated by (4‘1) has a suhsequenee {u^^)keN sueh that ^ 

in 00 and satisfies (3.12). 

Proof. Lemma [4.1| implies that 


- '^^^hiUh), sup Lh{u 

neNo 


n=0 


so that (u)()neN is bounded in and 


^ Je,hi}^h)i 


lim ||u”+i -u^ll = 0. 


(4.8) 
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In addition we infer from (3.5) and (3.13) that {y'f^)neN stnd {p'^)neN sire also bound ed in 
H^{^) and hence also in since dimVh < oc. In particular, we infer from (4.4) 

that the time steps Tn can be chosen to be bounded from below by a positive constant. As 
a result there exists a subsequence ^^nd (uh^yh^Ph) ^ ^ Vh ^ Vh such 

that 

^ Uh in y^^ y^^ p^^ Ph iii and a.e. in Q. 


In particular, y^ = Sh{uh) and ph satisfies (3.13). We finally deduce from (4.1) 


a 


[ 

Jn 


rik+l 

h 


D{-Vh 


u 


nk+l\ 


U 


■ (Vfe - <'= + ^))dx 


/ (a(v/i) - a(u 

Jn 


Uk + l') 
h 




^ 

^rik JQ 


'^fc+1 

h 


-<'=)• (v/, - + 


for all v/i G /C/i. Recalling (4.8) as well as > 7 we find that is a solution of (3.12) by 
passing to the limit k ^ 00 . □ 


5 Computational examples 

We use a preconditioned biconjugate gradient stabilized solver for the stationary forward 


problem (3.3) and the adjoint problem (3.13). To solve (4.1) we use the primal-dual active 


set method presented in cni, where the resulting system of linear equations is solved by 
applying the direct solver UMFPACK [ 2 T]. 

We set 

yobs = yh + ^n{x), ( 5 . 1 ) 

where n{x) is a random variable with the standard normal zero mean distribution, A G M 
and yh is the solution of 


/ a{uh)Vyh ■ Vxdx = / guxdo Mx e Vh 
Jo, Jan 

where Uh defines the objective curve. 

There is one regularisation parameter a. When the data is noisy we expect that a suit¬ 
able size of a is obtained by balancing the fidelity term with the regularisation term in 
the objective functional. The size of is determined by the need to obtain an accurate 
approximation of the regularised problem. We note that the thickness of the interfacial 
layer between bulk regions is proportional to s. In order to resolve this interfacial layer we 
need to choose h see [23] for details. Typically reasonable results are obtained with 
around 8 to 10 elements across the interface. Away from the interface h can be chosen 
larger and hence adaptivity in space can heavily speed up computations. In fact we use the 
finite element toolbox Alberta 2 . 0 , see 133, for adaptivity and we implemented the same 
mesh refinement strategy as in [5], i.e. a fine mesh is constructed for all variables 
and p'^ where 0 < < 1 for at least one index i G {1,..., r} and with a coarser mesh 

present in the bulk regions where = 0 or {u'^)i = 1 for all i G {1,... ,r}. In Figure 

we display a plot of the triangulation of 14 which illustrates the finer mesh within the 
interface. 
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In our computations we found it convenient to choose hmin — ^ as the minimal diameter, 
hmax = ^ as the maximal diameter of all elements and we set Tn = O.Ol/s. The stopping 
criteria we used to terminate the algorithm was the size of the residual to the first order 
optimality condition, i.e. ||(u^^^ — \x^)/Tn\\ < 1.0e“^. For each computation we state the 
number of iterations, L, required to reach this stopping criteria. 

In the case r = 2 we have U 2 — 1 — ui and the vector-valued Allen-Cahn inequality with 
two order parameters is reduced in the computations to a scalar Allen-Cahn inequality. 



Figure 1: A converged triangulation 


5.1 Results with r = 2 and d = 2 

In this section we see how our method compares with the one presented in m- In all the 
computations unless otherwise stated we set = (— 1 , 1 )^, Jj^^(F) := \\yr — 

s = 77 ^, ai = 3, a 2 = 0.5, a = 0.0001, A = 0.05 and 

lOTT 

{ —0.5 if X = —1 or 7 / = —1 
0.5 if X = 1 or ^ = 1. 

Figure [^displays the results we obtain when using the same initial curve (a circle of radius 
0 . 6 ) and objective curve (a ‘skinny’ ellitpse, x^/(0.07)^ + 2 /^/( 0 . 5 )^ = 1 ) that are used in 
Section 4.1 of m- In this simulation we set A = 0, as in m- The left hand plot in Figure 
[^displays the initial curve, the centre plot the objective curve and the right hand plot the 
computed solution u^. The number of iterations required to reach the stopping criteria was 
L = 4417. 

Figure takes the same form as Figure except that this time we compare our results 
with those displayed in Section 4.4 of m- The initial curve is again a circle of radius 0.6 
while the objective curve consists of two objects 

(a; +0.35)2 (y + 0.35)2 _ (x - 0.35)2 (y _ 0 . 35)2 ^ 

(0.25)2 + (0.3)2 (0.2)2 + (0.2)2 ’ 
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as in m we set A = 0. The number of iterations required to reach the stopping criteria 
was L = 11117. From this example we see that our phase field model successfully deals 
with topol^ical change. 

In Figure ^we plot the Residual - u^)/rn||, l\\Sh{uh) - VobsWo^ 

:= + i(l - |u,,p))d X and Je^hi'^h) versus iteration number, for the 

first 2000 iterations, for the computations displayed in Figures and From this figure 
we see that, in both computations, for the first 50 iterations there is a steep decrease in 
smd after that the decrease is much more gradual. We also see that the Residual 
decreases at a much slower rate than In Figure we display two intermediate 

results from the set-up in Figure the plots display after 150 iterations (left hand 
plot), after 500 iterations (centre plot) and after 4417 iterations, once the iteration has 
converged (right hand plot). From this figure we see that after 500 iterations the solution 
is approximating the shape of the objective curve reasonably well although the curve is not 
yet defined by a well defined interfacial region. 




Figure 2: Jj^^(F) := ||^r ~ 2 /o 6 s||^ 2 (^p initial curve (left hand plot), objective curve (centre 
plot), (right hand plot) 



Figure 3: Jfid{T) \\yr — 2 /c» 6 s||^ 2 (^p initial curve (left hand plot), objective curve (centre 
plot), (right hand plot) 


In Figure we follow the authors in Section 4.2 of m in seeing how noise effects the 
solution. We take the same initial and objective curves as in Figure and display the 
solutions obtained with A = 0.05 (left hand plot), A = 0.1 (centre plot) and A = 0.2 (right 
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Figure 4: Plot of Je.hiy^h) stud the Residual, versus the number of 

iterations: results in Figure(left plot), results in Figure(right plot) 



Figure 5: JfidiV) ||^r ~ 2 /o 6 s||^ 2 (^p after 150 iterations (left plot), after 500 
iterations (centre plot) after 4417 iterations (right plot) 
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hand plot). The number of iterations required to reach the stopping criteria were L — 4236, 
L — 4075 and L = 8941 respectively. 

In Figure we follow the authors in Section 4.5 of m in seeing how the value of the 
regularisation parameter a effects the solution. For the initial curve we take a circle of 
radius 0.7 and for the objective curve we take the ellipse j (0.5)^ + j (0.4)^ = 1. For the 
choice A = 0.05 we display the solutions obtained with a — 0.01 (top centre) a = 0.001 (top 
right) and a — 0.0001 (bottom left) a — 0.000025 (bottom centre) a — 0.0000025 (bottom 
right). The number of iterations required to reach the stopping criteria were L = 3918, 
L — 9183, L — 5550, L — 8441 and L — 21228 respectively. From this figure we see that 
(j — 0.001 and a — 0.0001 give the best approximations to the objective curve. 

In Figure we plot T/^^(u/i), Tj^^(u/i), Je,/i(u/i) and the Residual for the first 4000 it¬ 
erations, tor the computations displayed in Figure with a = 0.001, a = 0.0001 and 
a — 0.000025. From this figure we see that for a — 0.001 the initial decrease in Je^h{^h) is 
more gradual than for a = 0.0001 and a = 0.000025. 



Figure 6: JfidiV) := \\yr — 2/c»6s|obtained by taking A = 0.05 (left hand plot) 
A = 0.1 (centre plot) and A = 0.2 (right hand plot) 


In Figure]^ we show the effect that the size of |ai — a 2 | has on the solution u^. We display 
the objective curve in the left hand plot and in the subsequent plots we display a zoomed 
in image of the approximate solution, u^, at the end of the simulation obtained from 
decreasing value of \ai — a 2 |. We take a 2 = 0.5 in all plots and ai = 1, 3, 7 in the second, 
third and fourth plots respectively. We see that the approximation to the objective curve 
improves when |ai — a 2 | increases. The number of iterations required to reach the stopping 
criteria were L = 11891, L = 5550 and L = 17072 respectively. 

In Figure 10 we show the effect that the choice of O has on the solution u^. We com¬ 
pare results obtained by taking Jj^^(F) := \\yr — to results obtained by taking 

JfidiX) •= Wvt— yobs\\\ 2 (^QQy In thcsc simulations we set A = 0.02. We display the objective 
curve in the left hand plot and the approximate solution at the end of the simulation ob- 
tainedfrom JfidX) := \ \yr -yobs\\l 2 (^^~^ (centre plot) and JfidX) := \\yr-yobs\\l 2 (^QQ) (right 
plot). From this figure we see that the approximation to the objective curve obtained using 
Jfid{T) Wyr — yobsW^L^^dQ) offected more by the noise than the approximation that is ob- 
tained using J/jd(r) := ||yr-yobs|li 2 (Q)-Furthermore using J/jd(r) := ||yr-yobs|li 2 (Q) gives 
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Figure 7: Jfid{T) := \\yr — 2 /o 6 s||^ 2(^)5 objective curve (top left), obtained by taking 
a = 0.01 (top centre) a = 0.001 (top right) and a = 0.0001 (bottom left) a = 0.000025 
(bottom centre) a = 0.0000025 (bottom right) 





Figure 8: Plot of *7/^^(u/i), -^J^(u/i), Je^hi^h) ^md the Residual, versus the number of 
iterations for the results in Figure a = 0.001 (left plot) a = 0.0001 (centre plot), 
a = 0.000025 (right plot) 
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Figure 9: Jj^^(r) := ||^r ~ 2/o6s|objective curve (first plot), zoomed in plot of 
obtained by taking (ai,a 2 ) = (1,0.5) (second plot), (ai,a 2 ) = (3,0.5) (third plot) and 
(ai,a 2 ) = (7,0.5) (fourth plot) 


a better approximation to the objective curve than using Jj^^(r) := ||^r 

number of iterations required to reach the stopping criteria were L = 16151 and L = 18081 

respectively. 



Figure 10: Objective curve (left plot), obtained from Jj^^(F) := ||2/r ~ 2 /o 6 s||^ 2 (^) (centre 
plot) and Jfid{T) ||yr - 2/06^11^2(90) (right plot) 


In Figure [TT] we display results for three objective curves; we plot the objective curves in 
the upper row and the solution at the end of the simulation in the lower row. In these 
simulations we took a — 0.00001 and A = 0.005. The number of iterations required to 
reach the stopping criteria were L — 12720, L — 22296 and L — 36036 respectively. From 
this figure we see that our method results in good approximations of the objective curves. 


5.2 Results with r = 3 and d = 2 

In all the computations in this section we set = (—1,1)^, JfidiV) \\yr — 2 /o 6 s||^ 2 (^), 

s ai = 0.8, a 2 = 0.2, as = 0.3, a = 0.001, A = 0.0 and 

Stt 

0 if X = ±1 

—0.5 if 7/ = —1 

0.5 if ^ = 1. 
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Figure 11: J/id(r) := ||yr - yohs\\\ 2 ^^) 


, objective curves (upper plots), (lower plots) 


In Figure [T^ we display results for four objective curves, for each curve we took random t 
initial data for u^. We plot the objective curves in the upper row and the solution at 
the end of the simulation in the lower row. The number of iterations required to reach the 
stopping criteria were L — 10844, L — 33574, L — 31113 and L — 57373 respectively. 

5.3 Summary of the computational results 

The set-up of the computational examples presented in Figures Him and are taken 
from examples presented in m- The closeness of the approximated curve to the objective 
curve in the results that we present in Figures and is of a similar order to the results 
presented in m- In the case of Figure]^ the level set method used in m yields a better 
approximation to the skinny ellipse than our phase field model while in the case of Figure 
Hour results are a substantial improvement on the ones in jSI] as the level set method is 
unable to deal with the topological change required in this example whereas the phase field 
model successfully deals with it. 


6 Appendix 


Theorem 6.1. Let : X ^ M U {oo} be defined by 


Fein) := 



i/u e /C; 

otherwise. 


p 

Then F, where F is defined in (2.1). 
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Figure 12: Jj^^(r) := ||^r — 2 /o 6 s||^ 2(^)5 objective curves (upper plots), (lower plots) 


Proof. Let us first observe that for u G /C 


F,(u) = 


E/( 




+ 


i=l 


where : X := {v ^ L^(ri) | 0 < v{x) < 1 a.e. in 11} ^ M U {oo} is defined by 



if ^ G H^{n)nX; 
otherwise. 


It is well-known ([33], [I]) that ^ F with 


F{v) 


ly^d-i(^d*{v = 1} n , if G BV{n, {0,1}); 

8 

oo , otherwise. 


See (HldJllG] and the following development for the calculations leading to the factor tt/S. 
Let u G X and {uej^)keN C X an arbitrary sequence with lim^j^^oo — 0 and ^ u in 
Then C X and Ui in L^(ll), i = 1,..., r, so that 

r r r 

liminfF,^(ueJ = liminf > y^liminf> ^F{ui) = F{u) 

K rOO Ki rOO rOO 

i=l i=l i=l 


^ p ^ 

since F^ F. It remains to show that for every u G BV (H, {ei,..., e^}) H X there exists 
a sequence {ue}^)keN F fC with lim^j^^oo ^ such that ^ u in L^(H,M^) and 

limsupFefe(ueJ < F(u). (6.1) 

k^oo 
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We essentially follow the argument in [1]. Because of our particular choice of potential and 
the absence of volume constraints, the construction can be made more explicit allowing us 
at the same time to incorporate the condition that J2l=iUi{x) = 1 a.e. in which isn’t 
considered in [4j. 

Let u G BV{Q^{ei^... ,e^}) H X, say u = 'Y7i=iXEi^i- In view of Lemma 3.1 in |3] we 
can assume without loss of generality that the Ei are closed polygonal sets satisfying 
n = 0, i = 1,..., r. Lemma 3.3 in [3] implies that there exists rj > 0 such 
that the functions ^ M, 

[ dist(a;,9£;j), xeW^\Ei, 

hi{x) := < 

I —dist{x^dEi)^ X ^ Ei^ 


are Lipschitz-continuous on := {x G I \hi{x) \ < 7]} with \Vhi{x)\ = 1 a.e. in Let 
us introduce the function ^ 


(Pe{r) := < 


0 , 


1 

2 


(l + sin( 


T 

e 


1 , 


T < 0; 

|)), 0<r<e7r; 

T > evr. 


Furthermore, we define ^ ^ ^ by 


( l-iPeih) ,* = 1 ; 

ifeih) ■ ■ ■ (Pe{ti-i){l - (Pe{ti)) ,2 <i <r -1] 
(feih) ■ ■ ■ (pe{tr-l) ,i^r, 


where t = (ti,... It is not difficult to verify that 

Xe{t) = < 


ei , if < 0; 

Ci , if > 67r,.. ., ti-i > CTT, < 0; i = 2,..., r - 1; 

Cr , if > 67r,..., tr-1 > ctt; 

>r—l. 


0 < [Xe{t)]i <1, i = l,-■ ■,r \Dxe{t)\ <-&.e. in 

Xe{t) = i(l-sin(^-|))ei + i(l + sin(^-|))ej, 

if 0 < < CTT, tj < O^tk > CTT, A: = l,...,r—and i < j. 


(6.2) 


(6.3) 

(6.4) 


The above function is a particular example of the function constructed in Lemma 3.2 
in [3] . In addition we have 

r 

J2[xem = l, 

1=1 

As a consequence, the function Ue(x) := Xe{hi{x)^ ..., /i^_i(x)), x G belongs to JC and 
satisfies (see p. 79 in [3]) 

Ue ^ u in e ^ 0. 
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In order to analyze Fe(ue) we introduce as in [1] for i, j = 1,..., r the sets := £’i, 

:= {x e Ei\ hj{x) > enj = 1,... - 1}, i = 2,... ,r; 

^ij G 12 I 0 < hi{x) < 67r, hj{x) < 0, hk{x) > evr, A: ^ i, j} if i < j; 

Kl- := {x ^ Q\0 < hi{x) < ctt, 0 < hj{x) < en} if i < j. 

Then, 

r 

(6.5) 


1=1 


Kj 


Kj 


and 

Ue{x) = < 




X € 


1. /hjix) TT^S 1/ /hj{x) TT.. 


(6.6) 


X we have in view of (6.5) and (6.6) 


Abbreviating Fe(u, A) := f (^|T)up + ;^(1 — \u\‘^))d: 

J A 2 26 

FeAe) < ^ F,(u„ F,(u„ xy 

i<j i<j 

It is shown in |3j that limsup^^Q Fg(ue, = 0 for i, j = l,...,r, i < j. Furthermore, 
observing (6.6) and |V/i^(x)| = 1 a.e. in 12|^- we obtain 

|L>Ue(a;)|^ = ^cos2(—^--), 1 - |ue(a;)|2 =-cosy—^--), a; € %, 

SO that the coarea formula yields 

2 


(Ue,9^y = cosyhM _ l^dx = ^ ^ cosy^ - |) ^{{hi = t}n Ej)dt 


/ TT 

' cosys) n'^-\{hi = e(s + J)} n Ej)ds ^n'^-^{dEi n dEj n f2),e ^ 0. 

2 


4 


Hence, 


limsupF,(u,) < -J^'H'^-^dEindEjnn) = - n f2) = F(u), 

e^O 4 ^ 8 ^ 

1<J 1=1 


p 

where we note that dEiDdEj is counted twice in the second sum. In conclusion, E^ ^ E. □ 

Corollary 6.2. Suppose that (ue)e>o C E is a sequence such that {Ee{ue))e>ii is hounded. 
Then there exists a sequence e^ ^ ^ and u G BV{Tt^ {ei,..., e^}) H X such that ^ u 
in L^{n,W). 

Proof. Our assumption yields that (Te('^e,z))e>o is bounded for i = 1,..., r. It is well-known 
that this implies that there exists a sequence e^ ^ 0 and Ui G HH(12,{0,1}) such that 
Uej^^i ^ Ui in L^(12) and a.e. in 12, i = 1,..., r. Clearly, ^ u = (i/i,... ^Ur) in L^(12,M^), 
while it also follows that J2l=i — 1 i^^ ^ so that u G BV (12, {ei,..., Cr}) nX. □ 
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